{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum Likelihood Estimation of a multivariate Gaussian model\n",
    "\n",
    "The goal of this notebook is to use PyTorch to implement Gradient-based MLE for a multivariate Gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single 2D Gaussian component\n",
    "\n",
    "Let's generate a some data by sampling a multivariate Gaussian with an arbitrary covariance matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.49671415, -0.1382643 ])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng = np.random.RandomState(42)\n",
    "n_features = 2\n",
    "\n",
    "mean = rng.randn(n_features)\n",
    "mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.73912039, -0.50825603],\n",
       "       [-0.50825603,  0.10964792]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h = rng.randn(n_features, n_features)\n",
    "Cov = h @ h.T\n",
    "Cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.61547341,  0.        ],\n",
       "       [ 7.48827815,  3.01995037]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cholesky(np.linalg.inv(Cov))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHadJREFUeJzt3X+wXPV53/H3B3GxL/7BhXCHwAVZaqLiYERRfEtwVScG\nRCFNghScBly7xY0ZdSZhxmldmsswYzuOO9xY0zjp1DOtij1xYia2Y2OhRDjih+RJyyCHqwjMTxVC\nbaFrfsiA3GLJQT+e/nH3itXqnN2ze87Zc3b385q5o9295+75XiGe893n+3yfo4jAzMxGy0lVD8DM\nzPrPwd/MbAQ5+JuZjSAHfzOzEeTgb2Y2ghz8zcxGkIO/mdkIcvA3MxtBDv5mZiPo5KoHkObMM8+M\nZcuWVT0MM7OBsnPnzh9ExGSn42ob/JctW8bc3FzVwzAzGyiSvpflOKd9zMxGkIO/mdkIcvA3MxtB\nDv5mZiPIwd/MbAQ5+JuZjSAHfzOzEVTbOv+8Nu2aZ8PW3Xx//0HOmRjn5qvOZ92qqaqHZWZWC0MZ\n/DftmueWOx/l4KEjAMzvP8gtdz4K4AuAmRlDmvbZsHX3scC/6OChI2zYuruiEZmZ1ctQBv/v7z/Y\n1etmZqNm6IL/pl3znCQlfu+cifE+j8bMrJ6GKvgv5vqPRJzwvfGxJdx81fkVjMrMrH6GKvgn5foB\nlkjcdu1KL/aamTUMVfBPy+kfjXDgNzNrMlSlnudMjDOfcAHIm+v3ngEzGzZDNfO/+arzGR9bctxr\neXP9i+sI8/sPEryxZ2DTrvmcozUzq85QBf91q6a47dqVTE2MI2BqYjx3rt97BsxsGA1V2gcWLgBF\npmS8Z8DMhtFQzfzLkLZe4D0DZjbICgn+kq6WtFvSM5Jm2hz3fkkhabqI8/ZDGesIZmZVy532kbQE\n+BxwJbAXeEjS5oh4ouW4twEfBb6d95z9tJhCcrWPmQ2TInL+lwDPRMSzAJK+DKwFnmg57veA3wdu\nLuCcfVX0OoKZWdWKSPtMAc81Pd/beO0YST8LnBcRWwo4n5mZ5VT6gq+kk4A/AD6W4dj1kuYkze3b\nt6/soZmZjawigv88cF7T83Mbry16G3Ah8C1J3wUuBTYnLfpGxMaImI6I6cnJyQKGZmZmSYoI/g8B\nKyQtl3QKcD2wefGbEfHDiDgzIpZFxDJgB3BNRMwVcG4zM+tB7uAfEYeBm4CtwJPAVyPicUmfknRN\n3vc3M7PiFbLDNyLuBu5uee3jKce+r4hzmplZ74auvYN15i6lZubgP2IWu5QuNqtb7FIK+AJgNkLc\n22fEuEupmYFn/gOliHSNu5SaGXjmPzCKuqmMu5SaGTj4D4yi0jXuUmpm4LRPrTWneSLlmG7TNUV0\nKXW1kNngc/CvqdaqnDS9pGvydCl1tZDZcHDaJ6dNu+ZZPbuN5TNbWD27rbAbuyeleVqJheBb5Hl7\nGZerhcwGj2f+OZQ5C+6UzhEcSwX1c/btaiGz4eCZfw5lzoLT0jlTE+NMTYyfsAbQzXnzfFpxtZDZ\ncHDwz6HMWXC7qpw8581bMupqIbPh4OCfQ5mz4HWrprjt2pVMTYwjFmb8t127knWrpnKdN++nlXbj\nMrPB4Zx/Dpe9c5I7duw5LgVT5Cw4rSrn5qvOP6ESKOt5i/i04nsamw0+z/x7tGnXPF/fOX9c4Bfw\n/neXHxjzzL6dszcz8My/Z0npkwC2P9Wfew/3OvvO86nBzIaHg3+PBrXksYgdvmY2+Bz8e3TOxDjz\nCYF+ENInztmbmXP+PXLJo5kNMs/8e1S39ImbrZlZNwoJ/pKuBv4IWALcHhGzLd//98CNwGFgH/Ab\nEfG9Is5dpbqkT9xszcy6lTvtI2kJ8DngF4ELgA9IuqDlsF3AdERcBHwN+Eze89ob3GzNzLpVRM7/\nEuCZiHg2Il4HvgysbT4gIrZHxIHG0x3AuQWc1xoGtfLIzKpTRPCfAp5rer638VqajwDfLOC81pBW\nYRTQ13bPZjY4+lrtI+lDwDSwIeX76yXNSZrbt68/m6WGQVLl0aJe7/U7SMq6p4LZMCsi+M8D5zU9\nP7fx2nEkrQFuBa6JiL9PeqOI2BgR0xExPTk5WcDQRkNzu4ckw5z/L+rG9majpojg/xCwQtJySacA\n1wObmw+QtAr47ywE/pcKOKe1WLdqigdmLkcp3x/W/L8Xu816k7vUMyIOS7oJ2MpCqecXIuJxSZ8C\n5iJiMwtpnrcCfy4JYE9EXJP33Haiuu88Lno/ghe7zXpTSJ1/RNwN3N3y2sebHq8p4jzWWZ0bt5Wx\nH6HuFzuzunJ7hyFT55utlJGicZsNs964vcMQqsvO41ZlpGjq1mbDbFA4+BfI/XXaKytFU9eLnVmd\nOe1TkKJKDoe5Zt0pGrP6cPAvSBH57GGvWa/zeoTZqHHapyBF5LPbXUDKDJD9TFcNUorGaTwbZg7+\nBSkin11FzbrbQSfz34sNO6d9ClJEPjvtQlFmzbp3yCbz34sNOwf/ghSRz65iQbSOO2TrsOhdx78X\nsyI57VOgvPnsKmrWq9gh2y6XXpd0i3cO27Bz8K+Zfi+I9rsdRKfgXtWid6s6t8kwK4LTPiOu3+WX\nnXLpdUm3uCzVhp1n/pbr00a35ZCdgnud0i2DVJZq1i3P/K1nvWxK61TR5F3AZv3h4G8966UcMim4\nj50kDrx+mOUzW9iwdTfvf/eU0y1mJXPax3rWbX5+MUV08NARlkgciWBifIwfvX6YVw8cAhY+Pdyx\nYw8fvHQpn163sqvxtKagLnvnJNuf2ucdumYJPPO3nnWzKa05RQRwJOLYJ4BDR+K4YwO4Y8eerur7\nk1JQX9qxZ2j7JJnl5eBvPesmP5+WItp/8FDie0fjZ7JKev9W3qFr9ganfaxn3WxK66VUs5ufyXqs\nd+iaLXDwt1yylkOmlXCefuoY+w8cIlJ+JklSeWna+2d9T7NRU0jaR9LVknZLekbSTML33yTpK43v\nf1vSsiLOa4MjLUX0iV95Fx+8dClqOT4tfZRWXnrZOydPeP9WLhk1e0Pumb+kJcDngCuBvcBDkjZH\nxBNNh30EeDUiflrS9cDvA9flPbf1pqw+9e3et12KaN2qKabfcUamMaWtHWx/ah+3XbvS1T5mGSki\n6QN3F28gvQf4ZERc1Xh+C0BE3NZ0zNbGMQ9KOhl4AZiMNiefnp6Oubm5XGOzE7X21oGFGXHeWvqy\n3rfV8pktiSkiAf9n9pcKO4/ZoJK0MyKmOx1XRM5/Cniu6fle4OfSjomIw5J+CPwE8IMCzm9dKKtx\nWi/v2/xJ4bTxMSTYf+BQ21l6ndo/mA2yWpV6SlovaU7S3L59+6oezlAqq3FaLxu+mnP3+w8e4tXG\nwm+7mny3fzArRhHBfx44r+n5uY3XEo9ppH1OA15ufaOI2BgR0xExPTk5WcDQbNHiDVLS8mx5Z87d\n3oWsU11+Wk2+u22aFaOItM9DwApJy1kI8tcD/7LlmM3ADcCDwK8B29rl+61YSfn4ZkXMnLvtf5/l\nk8b8/oOsnt2WuECcNdj7JuxmyXIH/0YO/yZgK7AE+EJEPC7pU8BcRGwGPg/8qaRngFdYuEBYn7Sb\nZU+Mj/HJa96VKS/fLnh2exeyLHX5gmPH9HJHr7rcFcysjnJX+5TF1T7FSauQgYW0yQMzlyd+r8wK\nnk6fRgSJY2433larZ7clXmC6eQ+zQZO12qdWC75Wjnb5/Hbpl15aNmfVmrufGB/j9FPHjuXx0y5W\nRbR8cIsHMwf/kXDzVeefsIN2US8XhqKC57pVUzwwczmfve5i3vKmk48r85zqcgG5m2NdFmrm4D8S\n1q2a6qqFwqJ+BM9u2jV0uzDtslCzdA7+I+LT61by2esu7qpEsh/Bs1O7hjwlnS4LNUvnBV9rq4xS\nyeb3TPvX53YNZr3pZ3sHG2Ld1NRn0anKZ5Hz8mblcvC3vspyx61hz8t745nVgYO/9cViwGu3sUsw\n9MHQG8+sLhz8rXRZUj2jsvGqrK6qZt1y8LfSdUr1lJnmyZpi6VcqxhvPrC4c/K107QLbVImBNmuK\npZ+pGN+PwOrCdf5WurTAtkTi+/sPsmHr7sTe/XllbU9RZhuLVt54ZnXhmb+VLqndM8CRxh6Tsmba\n7VIsWfYalJGK6bb7qVlZHPytdK0B7yTpWOBf1M2iZ9b8fFqKZeLUsUr3GhS9d8KsFw7+1hfNAW/5\nzJbEY7LMtLvJz6fdYCaCUvYauH7fBomDv/VdnkXPtPz8Jzc/nhp4W1//d195uO05lki8/93dzc5d\nv2+DxsHf+q7bWz42S9sktv/gIfYfPHTsmObA2xp8O202OxLB13fOM/2OMzIHbtfv26BxtY/1Xa/d\nNjftmk+9L0GrdtU6SRU33fx8krSLiev3ra4887dK9LLouWHr7tTKnCRpgXfxvL/dIf2TNXAvXpSS\nxub6fasrz/xtYHQ7i24XeNetmkq9W1iWn2+WdlESuH7faitX8Jd0hqR7JT3d+PP0hGMulvSgpMcl\nfUfSdXnOaaMrLRiffupYTxun2qV/uqn2SbsoBV7stfrKO/OfAe6PiBXA/Y3nrQ4A/zoi3gVcDfyh\npImc57URlLY79hO/8q6e1hCa1x5gocqHLn5+UdpFqdMnC7Mq5bqTl6TdwPsi4nlJZwPfioi20yVJ\njwC/FhFPtzvOd/KyJL3W0pdZg5/UtXR8bIlvGWmV6NedvM6KiOcbj18AzuowqEuAU4C/S/n+emA9\nwNKlS3MOzQZNlgDdy0Jx2TX4i+/xu3/xOK8eWCg3fdPJXk6zeuv4L1TSfZIeS/ha23xcLHyESP0Y\n0fhk8KfAv4mIo0nHRMTGiJiOiOnJyckufxUbZIsBer7Ra2cxQBfR8K1fjdt+fOiNf9b7Dx4qbPxm\nZegY/CNiTURcmPB1F/BiI6gvBveXkt5D0tuBLcCtEbGjyF/AhkOZAbofPfT72RnUrAh5P5tuBm5o\nPL4BuKv1AEmnAN8A/iQivpbzfDakygzQaQuyRdbgp41zfv9Bls9sYfXstlI/BWzaNc/q2W19OZcN\nh7zBfxa4UtLTwJrGcyRNS7q9ccyvAz8PfFjSw42vi3Oe1wZMp+BUZoDuRw/9duMsOo3VqsyUmQ2v\nXME/Il6OiCsiYkUjPfRK4/W5iLix8fhLETEWERc3fbXfWmlDJUtwKjNA99pOohtltIzIyikn64Xb\nO1jpsjQ96/YmJ92WbpbdQ791/P28QYzvC2y9cPC30mUNTlkDdFrp5tz3XmH7U/sq66ffPP7Vs9v6\ndq9e3xfYeuFiZCtd0fn8tE8Sd+zYU5u8dzdprLyLtb4vsPXCwd9KV3RwatdLp1mVee+s6wxFLNb2\nY03Dho/TPla6om9aPnHq2LGdtJ20u2lL2bKksYq6CYzvC2zdcvC3vmgXnLpZvN20a57Xfnw483nV\n+JkqAmOW38uLtVYVB3+rVLd9dzZs3c2hoyfW0oyPncSPDx09IfUTjZ/pd/O3rL+XF2utKg7+Vqm0\ntMfHvvoIcOIFIG1GnBT4O/1MsyKavzVfPE6SONLSMTcpnZPnfsZmeXjB1yqVFpiPRCQufLarHErr\nn59lFp13o1Trwm1r4F+UVN7qxVqrgoO/VapdYE4Kvu0qh/JUFeXNvSddPJI4nWN14eBvlerUFqGb\nmXKeWXTevQhZLhJJFyL35bGqOOdvlVoMzB/76iOJqZKk4NuucqjXksc8ufdNu+YTc/ywcGvIoxGp\nC8hFlXqadcvB3yq3GOSqXPhstxehXRXQ4sw9KfBnuZVjGaWeZd6y0oaHg7/VQpaNYGUHtaRPDZ2q\ngNJy/UukTCmnoks9y75lpQ0PB3+rjU4bwaoIap3SMmkz9KMRiReS1otX0aWeTiNZVl7wtYFQVc/6\nTmmZrAvFaQu7QKGlnt4xbFl55m8Doaqg1iktk3Xm3u7i9cDM5YXNyr1j2LLyzN8GQj/uw5uk096B\nrOWl/bp4ub2zZeWZvw2EqtogZFmIzlJe2q8ZedEdVG14KVK2oVdteno65ubmqh6G1cgglzC2LljD\n8aWg3f5ug/x3YeWStDMipjsdl2vmL+kM4CvAMuC7wK9HxKspx74deALYFBE35TmvjaZB7lnfaR9B\nN5VMLue0IuSa+Uv6DPBKRMxKmgFOj4jfSTn2j4DJxvEdg79n/jYq0u73OzUxzgMzl+c+3kZLX2b+\nwFrgfY3HXwS+BZwQ/CW9GzgL+Cug46DMRkm3i8FZX29NDS37iXF2PPsqRyJYIvGBnzuPT69bmW/w\nNrDyBv+zIuL5xuMXWAjwx5F0EvCfgQ8Ba9q9maT1wHqApUuX5hya2WBIWww+bXyM1bPbTkgTZVk8\nTkoNNf/MkQi+tGMPgC8AI6pjqaek+yQ9lvC1tvm4WMgfJeWQfhO4OyL2djpXRGyMiOmImJ6cnMz8\nS5gNsqTyzLGTxI9eP5zY7TNLOWfWFtN/9u3njj3etGue1bPbWD6zhdWz29xZdMh1nPlHROpsXdKL\nks6OiOclnQ28lHDYe4D3SvpN4K3AKZJei4iZnkdtVpEyqmySFoMPvH74hJvUL97h7GgEp42P8eax\nk9h/4FDiOLLuH1hsSOdF5NGTN+2zGbgBmG38eVfrARHxwcXHkj4MTDvw2yAqM0C2VjItn9mSeNxi\nsN5/8BDjY0v47HUXJ547LTXUaokEuCfQKMq7w3cWuFLS0yzk82cBJE1Luj3v4MzqpJ/9hbJs/mp3\n7k43yVl06T84PbV6CNwTaJjlCv4R8XJEXBERKyJiTUS80nh9LiJuTDj+j13jb4Oqn/2FsgbvtHMn\ntZ1Y/VNnHJvpL5FY/VNn8Ld7ftj2E4J7Ag0vt3cwy6ifTdM63eEsy7k7bYpbPbut7aKwewINNzd2\nM8uo303T1q2a4mibwJ/33O0+seRtLW3155m/WUZVNE1L+7SR9U5hvby3dwqPBgd/sy70u79QWjfT\nImblVXVKtXpw8DersTI/bbj982hzS2czsyGStbGbF3zNzEaQ0z5mA6hdmwnf6MWycPA3GzDt2kwA\n7tFjmTj4mw2YTm0myuzR408Vw8PB32zA9NJmoogWFO78OVy84Gs2YNJaOpwzMd72e3n1s7Gdlc/B\n32zAtGszUWYLin42trPyOe1jNmCybM4qIy/fz8Z2Vj5v8jKzTFpz/lBcqwkrTtZNXp75m1kmbgcx\nXBz8zSyzfje2s/I4+JsNANfXW9Ec/M1qLq2+fu57r7D9qX2+IFhPcpV6SjpD0r2Snm78eXrKcUsl\n3SPpSUlPSFqW57xmoyStvv6OHXuY33+Q4I0LwqZd89UM0gZO3jr/GeD+iFgB3N94nuRPgA0R8TPA\nJcBLOc9rNjLS6uhb6/S84cq6kTftsxZ4X+PxF4FvAb/TfICkC4CTI+JegIh4Lec5zUZKWn19kjpv\nuPK6Rb3knfmfFRHPNx6/AJyVcMw/BPZLulPSLkkbJC1JOM7MEiTt2lXKsXXdcLW4btGcprr5zx9h\n1afuYfnMFlbPbnPKqs86zvwl3Qf8ZMK3bm1+EhEhKWnH2MnAe4FVwB7gK8CHgc8nnGs9sB5g6dKl\nnYZmNhKS6usve+ckX9853/X9d6uafSetWxw6Grx64BDgJnFV6Bj8I2JN2vckvSjp7Ih4XtLZJOfy\n9wIPR8SzjZ/ZBFxKQvCPiI3ARljY4ZvtVzAbfkn19dPvOKOrQF5mV85OF5Us6agiW09bZ3lz/puB\nG4DZxp93JRzzEDAhaTIi9gGXA+7bYJZTtxuu2nXlzBNws1xUsq5b1HnNYtjkzfnPAldKehpY03iO\npGlJtwNExBHgPwD3S3qUhXTl/8h5XjPrUlldObO0ek5at0hS1zWLYZRr5h8RLwNXJLw+B9zY9Pxe\n4KI85zKzfMrqypnlotK6bnHa+Bg/ev0wh468kd0tqvW0ZeMdvmYj4uarzk/sypk34Ga9qLSmqVz6\nWS0Hf7MRUVZXzl4vKm4SVy0Hf7MRUkbAdavnweTgb2bHyZqOcdpmsDn4m9kxWfcCdDrOF4b6c/A3\ns2Oy7gXoVN5Z1mYyK07eOn8zGyJZ9wK0Oy5L3b9Vz8HfzI5Jq/lvfb3dcWVtJrNiOfib2TFJO3GT\nyjbbHZf1AmLVcs7fzI7JWrbZ6bgyNpNZsRRRz+aZ09PTMTfn/m9mg8jVPtWRtDMipjsd55m/mRWu\nqt27vuhk5+BvZrXTSxAv834Fw8gLvmZWK0m3fLzlzkc73ubRJabdcfA3s1rpNYi7xLQ7Dv5mViu9\nBnGXmHbHwd/MaqXXIJ51j4ItcPA3s1rpNYivWzXFbdeuZGpiHAFTE+Pcdu1KL/amcLWPmdVKnvsD\n+AYx2Tn4m1ntOIiXL1faR9IZku6V9HTjz9NTjvuMpMclPSnpv0hSnvOamVk+eXP+M8D9EbECuL/x\n/DiS/gmwGrgIuBD4x8Av5DyvmZnlkDf4rwW+2Hj8RWBdwjEBvBk4BXgTMAa8mPO8ZmaWQ97gf1ZE\nPN94/AJwVusBEfEgsB14vvG1NSKezHleMzPLoeOCr6T7gJ9M+NatzU8iIiSd0CJU0k8DPwOc23jp\nXknvjYj/mXDsemA9wNKlSzuP3szMetIx+EfEmrTvSXpR0tkR8byks4GXEg77VWBHRLzW+JlvAu8B\nTgj+EbER2AgLLZ2z/QpmZtatvGmfzcANjcc3AHclHLMH+AVJJ0saY2Gx12kfM7MK5Q3+s8CVkp4G\n1jSeI2la0u2NY74G/B3wKPAI8EhE/EXO85qZWQ65NnlFxMvAFQmvzwE3Nh4fAf5tnvOYmY2Cft6M\nxjt8zcxqoN83o3FjNzOzGuj3zWgc/M3MaqDfN6Nx8Dczq4F+34zGwd/MrAb6fTMaL/iamdVAnvsY\n9MLB38ysJvp5HwOnfczMRpCDv5nZCHLwNzMbQQ7+ZmYjyMHfzGwEOfibmY0gB38zsxGkiHreMEvS\nPuBHwA+qHksbZ1Lv8YHHWIS6jw/qP8a6jw+GZ4zviIjJTm9U2+APIGkuIqarHkeauo8PPMYi1H18\nUP8x1n18MHpjdNrHzGwEOfibmY2gugf/jVUPoIO6jw88xiLUfXxQ/zHWfXwwYmOsdc7fzMzKUfeZ\nv5mZlaDWwV/S70n6jqSHJd0j6Zyqx9RK0gZJTzXG+Q1JE1WPqZWkfyHpcUlHJdWmmkHS1ZJ2S3pG\n0kzV42kl6QuSXpL0WNVjSSLpPEnbJT3R+O/70arH1ErSmyX9jaRHGmP83arHlETSEkm7JP1l1WNJ\nIum7kh5txMK5It6z1sEf2BARF0XExcBfAh+vekAJ7gUujIiLgP8N3FLxeJI8BlwL/HXVA1kkaQnw\nOeAXgQuAD0i6oNpRneCPgaurHkQbh4GPRcQFwKXAb9Xw7/Dvgcsj4h8BFwNXS7q04jEl+SjwZNWD\n6OCyiLh4JEo9I+L/Nj19C1C7BYqIuCciDjee7gDOrXI8SSLiyYjYXfU4WlwCPBMRz0bE68CXgbUV\nj+k4EfHXwCtVjyNNRDwfEX/bePz/WAhe/bkTSEax4LXG07HGV63+P5Z0LvBLwO1Vj6Wfah38AST9\nJ0nPAR+knjP/Zr8BfLPqQQyIKeC5pud7qVngGiSSlgGrgG9XO5ITNVIqDwMvAfdGRN3G+IfAfwSO\nVj2QNgK4R9JOSeuLeMPKg7+k+yQ9lvC1FiAibo2I84A7gJvqOMbGMbey8DH8jrqO0YaTpLcCXwd+\nu+XTci1ExJFG6vZc4BJJF1Y9pkWSfhl4KSJ2Vj2WDv5pRPwsC2nS35L083nfsPJ7+EbEmoyH3gHc\nDXyixOEk6jRGSR8Gfhm4Iiqqne3i77Eu5oHzmp6f23jNuiBpjIXAf0dE3Fn1eNqJiP2StrOwjlKX\nRfTVwDWS/jnwZuDtkr4UER+qeFzHiYj5xp8vSfoGC2nTXGt4lc/825G0ounpWuCpqsaSRtLVLHxk\nvCYiDlQ9ngHyELBC0nJJpwDXA5srHtNAkSTg88CTEfEHVY8niaTJxQo4SePAldTo/+OIuCUizo2I\nZSz8G9xWt8Av6S2S3rb4GPhnFHDxrHXwB2YbqYvvsPAL166UDfivwNuAextlWP+t6gG1kvSrkvYC\n7wG2SNpa9Zgai+Q3AVtZWKj8akQ8Xu2ojifpz4AHgfMl7ZX0karH1GI18K+Ayxv/9h5uzGDr5Gxg\ne+P/4YdYyPnXspyyxs4C/pekR4C/AbZExF/lfVPv8DUzG0F1n/mbmVkJHPzNzEaQg7+Z2Qhy8Dcz\nG0EO/mZmI8jB38xsBDn4m5mNIAd/M7MR9P8BtE8yyzbe98gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3b1979f128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_samples = 100\n",
    "data = rng.multivariate_normal(mean, Cov, size=n_samples)\n",
    "plt.scatter(data[:, 0], data[:, 1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the MLE estimate from this data using the closed-form formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.70263098, -0.17305392])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_mle = data.mean(axis=0)\n",
    "mu_mle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.0476404 , -0.38824024],\n",
       "       [-0.38824024,  0.08873728]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Cov_mle = (data - mu_mle).T @ (data - mu_mle) / n_samples\n",
    "Cov_mle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.69266713,  0.        ],\n",
       "       [ 7.40569797,  3.356966  ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cholesky(np.linalg.inv(Cov_mle))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parametrisation of a positive definite matrix\n",
    "\n",
    "\n",
    "Let's parametrize the precision matrix `P` (inverse of a covariance matrix `C`) as follows:\n",
    "\n",
    "- `P` has Cholesky decomposition `H`\n",
    "- `H` is a lower triangular with a positive diagonal\n",
    "- the log of the diagonal entry is stored in a vector of parameters named `d`\n",
    "- the off diagonal elements of `H` are stored in the matrix of parameters named `W`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 7.3891 -1.4755\n",
       "-1.4755  7.6837\n",
       "[torch.FloatTensor of size 2x2]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import torch\n",
    "from torch.autograd import Variable\n",
    "\n",
    "mu = Variable(torch.zeros(n_features), requires_grad=True)\n",
    "d = Variable(torch.ones(n_features), requires_grad=True)\n",
    "W = Variable(torch.randn(n_features, n_features), requires_grad=True)\n",
    "H = torch.diag(torch.exp(d)) + torch.tril(W, -1)\n",
    "P = H @ H.transpose(1, 0)\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check that H is the actual Cholesky decomposition of P:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.71828175,  0.        ],\n",
       "       [-0.54281157,  2.71828175]], dtype=float32)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H.data.numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.71828175,  0.        ],\n",
       "       [-0.54281157,  2.71828175]], dtype=float32)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cholesky(P.data.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`P` is positive semi-definite by construction (product of a matrix `H` by its transposed).\n",
    "\n",
    "Because of we take the `exp` of `d` to build the diagonal elements of `H`, the determinant of `H` and therefore `P` is stricly positive.\n",
    "\n",
    "`P` is therefore is positive definite, whatever the values the parameters in `d` and `W`. Because the Cholesky decomposition exists for any symmetric positive-definite  matrix and is unique and `exp` is a bijection from $\\mathbb{R}$ to $\\mathbb{R}^+$, this parametrization of the manifold of positive definite matrices is bijective."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "54.598145"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.det(P.data.numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "54.598144564972472"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.det(H.data.numpy()) ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The determinant of `P` is cheap to compute from the `d` parameters directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "54.59814456497247"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "torch.prod(torch.exp(d.data) ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the above function to define the log-likelihood of a Gaussian model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from math import log\n",
    "\n",
    "\n",
    "def loglik(X, mu, d, W):\n",
    "    \"\"\"Compute the average log-likelihood of samples\n",
    "    \n",
    "    X shape: (n_samples, n_features)\n",
    "        data points\n",
    "        \n",
    "    mu: shape: (n_features,)\n",
    "        parameters of the mean of the Gaussian.\n",
    "    \n",
    "    d: shape: (n_features,)\n",
    "        parameters of the diagonal of the Cholesky factor of the\n",
    "        precision matrix of the Gaussian.\n",
    "        \n",
    "    W: shape: (n_features, n_features)\n",
    "        parameters of the off-diagonal of the Cholesky factor of the\n",
    "        precision matrix of the Gaussian. The upper-diagonal elements\n",
    "        are ignored.\n",
    "    \"\"\"\n",
    "    H = torch.diag(torch.exp(d)) + torch.tril(W, -1)\n",
    "    P = H @ H.transpose(1, 0)\n",
    "    diff = X - mu\n",
    "    quad_form = torch.sum(diff * (diff @ P), dim=1)\n",
    "    return (-0.5 * log(2 * np.pi) + torch.sum(d) - 0.5 * quad_form)\n",
    "\n",
    "\n",
    "def nll(X, mu, d, W):\n",
    "    \"\"\"Average negative log likelihood loss to minimize\"\"\"\n",
    "    return -torch.mean(loglik(X, mu, d, W))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 9.5162\n",
       "[torch.FloatTensor of size 1]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = Variable(torch.FloatTensor(data))\n",
    "loss = nll(X, mu, d, W)\n",
    "loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "loss.backward([torch.ones(1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       "-5.4471\n",
       " 2.3664\n",
       "[torch.FloatTensor of size 2]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu.grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 18.5303\n",
       " -0.1230\n",
       "[torch.FloatTensor of size 2]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d.grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 0.0000  0.0000\n",
       "-1.4503  0.0000\n",
       "[torch.FloatTensor of size 2x2]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W.grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 0.1816\n",
       "[torch.FloatTensor of size 1]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_mle = torch.FloatTensor(np.linalg.cholesky(np.linalg.inv(Cov_mle)))\n",
    "d_mle = Variable(torch.log(torch.diag(H_mle)))\n",
    "W_mle = Variable(torch.tril(H_mle, -1))\n",
    "nll(X, Variable(torch.FloatTensor(mu_mle)), d_mle, W_mle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 9.516191482543945\n",
      "100 0.36440008878707886\n",
      "200 0.19654464721679688\n",
      "300 0.1821846216917038\n",
      "400 0.18160586059093475\n",
      "500 0.18159542977809906\n",
      "600 0.18159565329551697\n",
      "700 0.18161100149154663\n",
      "800 0.18159569799900055\n",
      "900 0.18159551918506622\n",
      "1000 0.18159547448158264\n",
      "1100 0.18226905167102814\n",
      "1200 0.18159545958042145\n",
      "1300 0.1815958023071289\n",
      "1400 0.18162307143211365\n",
      "1500 0.18159577250480652\n",
      "1600 0.1815953552722931\n",
      "1700 0.18165569007396698\n",
      "1800 0.18159538507461548\n",
      "1900 0.18159548938274384\n"
     ]
    }
   ],
   "source": [
    "learning_rate = 1e-1\n",
    "optimizer = torch.optim.Adam([mu, d, W], lr=learning_rate)\n",
    "for t in range(2000):\n",
    "    # Compute and print loss.\n",
    "    loss = nll(X, mu, d, W)\n",
    "    if t % 100 == 0:\n",
    "        print(t, loss.data[0])\n",
    "    optimizer.zero_grad()\n",
    "    loss.backward()\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.70263098, -0.17305392])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_mle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 0.7028\n",
       "-0.1725\n",
       "[torch.FloatTensor of size 2]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 0.5263\n",
       " 1.2110\n",
       "[torch.FloatTensor of size 2]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_mle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 0.5263\n",
       " 1.2110\n",
       "[torch.FloatTensor of size 2]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       "-0.0000  0.0000\n",
       " 7.4055 -0.0000\n",
       "[torch.FloatTensor of size 2x2]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W * Variable(torch.tril(torch.ones(n_features, n_features), -1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 0.0000  0.0000\n",
       " 7.4057  0.0000\n",
       "[torch.FloatTensor of size 2x2]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W_mle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixture of 2D Gaussian components\n",
    "\n",
    "Let's generate some ground truth Gaussian Mixture Model with 3 components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([ 1.49014246, -0.4147929 ]),\n",
       " array([ 1.94306561,  4.56908957]),\n",
       " array([-0.70246012, -0.70241087])]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng = np.random.RandomState(42)\n",
    "n_features = 2\n",
    "n_components = 3\n",
    "\n",
    "true_means = [rng.randn(n_features) * 3 for _ in range(n_components)]\n",
    "true_means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "true_covariances = []\n",
    "for _ in range(n_components):\n",
    "    h = rng.randn(n_features, n_features)\n",
    "    true_covariances.append(h @ h.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's generate some data from the ground truth model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+cVPV5L/DPM7uz7LDKLggRWEAwMdooEHBjvGpMDb0Q\nuyqUJmiT+Eqa1yvctmlCba4VtcWNNxGU3hjSprcvr9rmldgqEoLGbQqp2tjoVbuIrjFKE7EIy6JY\n2EWXgZ3dee4fZ8/uzJnv+TVzdn6c/bz/iO7ZM+ecIfjMd57v832+oqogIqL4SFT6AYiIKFoM7ERE\nMcPATkQUMwzsREQxw8BORBQzDOxERDHDwE5EFDMM7EREMcPATkQUM/WVuOn06dN1/vz5lbg1EVHN\n2r179zuqOsPvvIoE9vnz56Orq6sStyYiqlkisj/IeUzFEBHFDAM7EVHMMLATEcUMAzsRUcwwsBMR\nxQwDOxFRzDCwE8VZ91bg7guAjhbrn91bK/1EVAYVqWMnojLo3gr8+KtAJm393H/A+hkAFq2p3HPR\nuOOInSiuHr99LKjbMmnrOMUaAztRXPUfDHecYoOBnSiumueEO06xwcBOFFfLNgDJVP6xZMo6TrHG\nyVOiuOneauXR+w8CqalAfQpIH7NG6ss2uE+c5r7O71yqagzsRHHirIRJH7VG6avv8Q7S3VuBR74M\nDA9aP/cfsH4GGNxrEFMxRHFSbCXMT24aC+q24UHrONUcBnaiOCm2EiZ9NNxxqmoM7ERxUslKGK5y\nrRoM7ERxUmwlTGpauONOdm6//wAAHVvl+tifMthXAAM7UZwsWgNc/R2geS4Asf559Xf8J0CvvBNI\nJPOPJZLW8SDccvtd9xcGewb3cceqGKK4WbQmfCWLfX6x5Y6uOXzN/9GeyGWlzbhiYCciSzEfCLbm\nOSMj8wDY0mDcMRVDRKUz5fYh5nPZ0mDcMbATUelMuf22L7KlQYUwFUNE0TClcuZdzDYFFcDATkTj\np5S8PRUtklSMiNwgIq+IyC9E5B9FpDGK6xIRUXglB3YRaQXwVQBtqnoBgDoA15V6XSIiKk5Uk6f1\nAFIiUg9gMoBDEV2XiKoNWwdUvZIDu6r2APhLAG8C6AXQr6q7Sr0uEZVgvIKvW+sABveqEkUqZiqA\nlQAWAJgNoElEPmc4b62IdIlI15EjR0q9LRG5MQXf7V8C7lxQegDmBtk1IYpUzG8BeENVj6hqBsB2\nAJc4T1LVe1S1TVXbZsyYEcFticjIFHwBqwVvqaNrbpBdE6II7G8CuFhEJouIAFgG4NUIrktExfAK\nsqWOrrlBdk2IIsf+HIBtAF4A8PLINe8p9bpEVCS/IFvK6JobZNeESKpiVPU2VT1PVS9Q1etV9VQU\n1yWiIhj7tuQIOro2TcAW2xaYyoorT4nipHtrTo5dUNA2N+jo2rkptl39AnA1aQ1gEzCiuMirhgEA\ntTbLSE1D6NE1q19qGkfsRHFhCsbZDNDQBNz0RrhrsfqlpnHEThQXUQbjoNUvXIValRjYieIiylLE\nINUvXIVatRjYieIiylLEINUvzMNXLebYieLCb0Pq7q3AT26yVqAC1qTqlXfmB2u7qibIxhjMw1ct\nBnaiOHErRezeCuz4I2sy1ZY+Cjzy5bHX+ZU4OrltYM1VqBXHVAzRRPD47flB3TY8OJY6CZta4SrU\nqsXATjQReKVH7FF32NQKV6FWLaZiiCYCt7QJAECsNEwxqRWuQq1KHLETTQTLNlirUI3UmlQdHCj8\nFVMrNYkjdqKJwB5Vb/+S+fd2pUwuZ9VMmIoZqiiO2IkmikVrRvLhATU05Qd1LkaqGQzsRBOJX0vf\nXLmTplyMVFMY2IkmElMlS2qa+dzcSVMuRqopDOxEceds1AVYI/fmOWOB2Tmx6pw05ZZ4NYWTp0Rx\nZlpN+siXAdWxBUvpo0BdgzVyTx8zT4wu25B/HYAVM1WMgZ0ozky58eHBwvOGB737tvv1oaGqwsBO\nFGdhcuD9B6xUjVvg5mKkmsEcO1Gchc2Bs5wxFhjYieIsTHmjE8sZaxZTMUS1wG3Vp99qUGduHBru\nvv0HgK9PA3TYKo1kXr0miGrI/6NNFxFpAXAvgAtg/c35oqr+P7fz29ratKurq+T7Ek0IzsoWwBqF\nL/4M8NI/FB736rB49wUezcAC8Ls+jSsR2a2qbX7nRZWK2QLgn1X1PACLAbwa0XWJyG3V5+6/918N\n6qxhP2d58akZ0/WpKpUc2EWkGcDlAO4DAFUdVNW+Uq9LRCPcKlt02Pt8U3+XrvuBOReNrTxNTbNq\n2KN4HqoaUYzYFwA4AuDvRGSPiNwrIk3Ok0RkrYh0iUjXkSNHIrgt0QThVtkidd7nm0b6UOCNp6xc\neUefVbe+8rtjgd7tmkGeh6pGFIG9HsBSAP9HVZcAGACw3nmSqt6jqm2q2jZjxowIbks0QbhtQXfh\nF7y3pnMdWWt+OmXRGuCGX1iBXrPez8LVpjUhisB+EMBBVX1u5OdtsAI9EUXBbQu6q77lvTWd18ja\nLeinprq/hlvf1YySyx1V9bCIHBCRc1V1L4BlAH5Z+qMR0Si3VZ9eq0GXbQC2r4WxxNEU9Lu3AoPv\nFR5PJIFVf8OAXkOiqor5CoAHRKQbwIcB3BHRdYmoWIvWAAsuN//unOWFxx6/3dxHZtLpDOo1JpIF\nSqr6IgDf2koiKrOj+8zHf7Wr8JhbeiZ9LLrnobLgylOiuDCtQg2zQUbzHPPiJa+8O1Ul9oohigO3\nPUndgrIpx75sg7mm/dS7+c3AnIue2Cis6jCwE9Wy7q3AnQuA7V8yr0IFvEsicy1aU7iTEmBtyGGX\nR3JT65rAwE5Uqx77Uyugp4+6n5M+5l0Smat7K5AZMF/HTt1wU+uawMBOVIu6t1rtAfykpgbf9cgr\nONupG25qXRMY2Ilq0eO3w7cFb12DlR8PmjbxCs526oabWtcEBnaiWhRkhJxIjm1YbfNKm7gF59S0\nsVG+W3sDthmoKgzsRLUoyAjZL1/u5Ba0r7xzrBJm+1qgPmUFe7s7ZH3KOs4KmarBwE5Ui0rZ8s6t\nBNKtJw2QXwmTPgoMpYG2L1r/TB8FK2SqCxcoEdUiOzXyo//h3pFREv7dGk3XdU6u3n2B+0Yfzp7w\ndqqHLQgqioGdqBY4V5Wesxx45Ufegdvtd2FbBBS70QdVDAM7UbVz7nnafwDous/7NckmYPI09xYB\nd18QrAQScG81IHXm4M4KmYpjjp2o2hl3QvKRzZj3N00krda8uSWQ278EfHO2e2682I0+qGIY2Imq\nXTGpjeFBq4OjczJ00unm1ryZAWDHH5mDe7EbfVDFiKrPIodx0NbWpl1dXWW/L1FNunOBd9sAT5Kf\nbulo9j69ea61TR5VJRHZraq+LdKZYyeKNTvdshbY8wMAAs8Vq5z4jAUGdqJqF8lGFwq88TP/0zjx\nGQvMsRNVO69gK3XR3SeR5MRnTDCwE1U7tw0wEkngd/52ZPKyRKlp3LA6RhjYiaqNc4ciAFj53ZH+\nLCOSTVaFy/a1wOCAOfAHkUwBq/8vcNMb1s/cGSkWGNiJqonbDkWAFXw7+q1AjOxYj5b0UUB1rDFX\n81xgwcetf/ckwOLPWKN07owUKwzsRNUkyA5FpnOyGaChCejos8oVP/8osPoen5upVese9L5UMyIL\n7CJSJyJ7ROSxqK5JNOEE2aEo6C5Gi9b459/7D1qjclPLAK97UVWLcsS+DsCrEV6PaOIJskNRmF2M\n/Nr7pqaOpXrCPA9VtUgCu4jMAdAO4N4orkc0YQXZoSjMLkZ2O4Dcidfc1wDufWjY96VmRTVi/zaA\nPwMQsvkzEeVx68uSW4YY5BznNW96w5p0db7Ga/ET+77UrJJ7xYjIVQB+W1X/SER+E8D/VNWrDOet\nBbAWAObNm3fh/v37S7ovEUXg7gvM+XX2jKlKQXvFRDFivxTANSLynwAeBPAJEfmB8yRVvUdV21S1\nbcaMGRHclmiCcNa1u5UgBj0vFzenjqWSA7uq3qyqc1R1PoDrADyhqp8r+cmIKHh9ebF16GHTOlQT\n2ASMqJp51ZfnBt+f3GQ+7yc3jV3Hbcck0z6nVNMiDeyq+q8A/jXKaxJNaEFq1ru3uvdrTx8FHvny\n2OYa9kj+zWetxUlBt8ejmsKVp0TVLEjNut/qUOeOSZk00HU/2wfEGAM7UTULMrlZ1OpQRzUc2wfE\nCgM7UTULMrkZ1epQtg+IDU6eElU7v8nNZRusVIppBWkiCSALZIf978P2AbHBETtRLTHVqueN6jG2\nq1LzXGvzjLpG/+uydj1WOGInqhV2rbo9MjdWuMzNr3Dp3gpkBjwuKqyKiSEGdqJa4VbT3nXf2M8F\nwd6lHS/AtgExxsBOVI26t1qLi+z69GSTz8g7h13O6Kx8cWLqJbYY2ImqTffW/EVFQPCgPsonqKem\nMfUSY5w8Jao2j99euKgoSskUcOWd43d9qjgGdqJqM5715GzyNSEwFUNUbt1bvZtyNc/xnvQMIpG0\nNri2JVPhArrfM1JV44idqJyCtNddtgGoayjtPiIj2+EV0Yq32BbAVDUY2InKyasNr23RGmDld/P3\nKU02ARLiP9fhQaChCejos0oaw4y2gzwjVTWmYojKKUgbXsDcRqCjJZp7Ffs69pKpGRyxE5VTkDa8\nTnYbAb8SxjDXLOZ17CVTMxjYicop7B6jefnuIu5VDO6DWvMY2InKKeweo6Z8t615bn4ePlcpC5C4\nD2rNY46dqMx2DF+Kzae+g0Mn05jdmMKNw+dildvJrnltsSZFnY3BgGgWIHEf1JrGETtRGe3Y04Ob\nt7+Mnr40FEBPXxo3b38ZO/b05J/ol1e3890cXZMBR+xEZbR5516kM/mbXqQzw9i8cy9WLWm1DphG\n4bmc+W6OrsmBI3aiMjrUZw7Wecf98uockZMPjtiJymh2Swo9huA+uyWnCsUvr07ko+QRu4jMFZEn\nReSXIvKKiKyL4sGI4ujGFecilazLO5ZK1uHGFeeOHSiljty0dR5NOFGkYoYAfE1VPwTgYgBfFpEP\nRXBdotrlEmBXLWnFxtUL0dqSggBobUlh4+qFY/l1oPg6cvZ4oRElp2JUtRdA78i/vysirwJoBfDL\nUq9NVJPc9iYFgEVrsGpJa34gd7Lz52G7K3r1eGFOfkKJNMcuIvMBLAHwnOF3awGsBYB58+ZFeVui\n6hJFgC2m0iVMjxe25Y21yKpiROQ0AD8E8Ceqetz5e1W9R1XbVLVtxowZUd2WqPpUqolW0Nw8Uzax\nF0lgF5EkrKD+gKpuj+KaRDUrzORnlJOdHrn5zn2dWL5tORZ9bxGWd92OzgbJP49teWMliqoYAXAf\ngFdV9VulPxJRjQs6+Rn1yNllFWrnaU3oeKYDvQO9UCh66wQd06ehs2ly/uvZljc2osixXwrgegAv\ni8iLI8duUdV/iuDaRLUnyORn91bgR38AaP4q1JInOw25+S3bluPk8Mm8YycTCWyZ2oL2gRNjB9mW\nNzaiqIr5OQDxPZFoIvGa/LRH6s6gbot45Hx44LD5eH1OPT3b8sYKV54SlZtXywCg9JGzo+Jl5pkt\n6M30F5w2MwtYKRtWxcQNAztRuXmNyEOMnDv3dWLLC1tweOAwZjbNxLql69D+3kBBDf06fRcd08/A\nSc2MvlazSRw7fh12rPysd0091SQGdqJya55j3hFJ6gI3+Orc14mOZzpGc+e9A73oeKYD6E+j3fFt\noP14HzDpdGyceib6Bt+GZlpw6sgKvHf8fNy8/WUAYHCPGXZ3JCo3t6qZ3/nbwOmQjc9+q3BCdPgk\ntkwy5+3bjxyEvnkr3nttEwZeX4+h40sAjLUMrjZ55ZnblqNzX2elH6mmcMRONE527OnB5p17cagv\njdktKdy44lxrZFxsy4Cc6/YNvg0xlCz01tcVHgSA5jk49FaAlsEhGdNBZ7cXfT37msZvI0DJ154o\nGNiJxoG9U5K9qYa9UxKAseBe5GTl5p17gZkC0+5KCRFYRWo5vxvJ209/4hWcaPoxJNk3mo4ZOr4k\nv2VwCOMVgLe8sMX8beSFLQzsATEVQxTSjj09uHTTE1iwvhOXbnqicFs7eO+UVCprhG3eMi+r6vid\nAIs/g87TmjA8bSsSDX0QARINfWictR2Tp76U3zI4BK8AXArX8kyX41SIgZ0ohKB7lgbaKalIs1tS\n0EyL8Xezhpw5dgV+tQtbXtiCjJ7K+40kMpg29/GiJ07HKwDPbJoZ6jgVYmAnCiHoSNwtvVFs2iPX\njSvOhR69EppN5h1vzGax7lhf4Qv6D7oG2+OZI0U/x3gF4HVL16GxrjHvWGNdI9Yt5R4+QTGwE4UQ\ndCQeaKekkOxKkQ3dV2La3MeRGJwPVQEUECSwclDyWwTYmuegeVKz8ZqNdY2+1SduFSrjFYDbz25H\nxyUdmNU0CwLBrKZZ6Likg/n1EDh5ShRCoD1LMVYXbqyKKYJzorI/8zbQ+PZoLw9FFo+kGrBkSotV\nt25LptC55Hfw3puPGa+bHk4jPWC9H9PkZ5AJ0qirYuxrM5AXT1TNkzDjqa2tTbu6usp+X6JSOatd\nAGskXrC9XcSWb1uO3oFe3/NmJZux662+vDLK5f9xb6DXjl6jaRZ2fWqX531zz8k1HuWPNEZEdqtq\nm995HLEThTSpPjEa2KdOTuK2q88vKqi71rkbBJ2QPJw5Dtzwi/xje74R6rly79UbYoKU9efVgzl2\nooDs0XpfeqznyslMtqRr5VbX/PxHf4MTd55n3HQj6ITklGTh7mRhJzPt83fs6YFmzLl50zXHq/yR\nwmNgpwktSE26raTa9JGdkjo3z8by+y/AX7z0SSTmfRP1U/YAAK5J/By3yz2YnO6FadMN00SlM4uq\n2SSOHlhW8B5Mr3WTO/m5eedenHx7RUH1DbJJ4wQp68+rB1MxNGH5rg51CFoR40yxfPtDv8JHXr4N\nnQ3WzkUnE9aUZ6KhD42zH4Ke+WO0/dcxTD4xmH/hTBonfrIBkxetGU1lbHz2W6ONvIbeOw/1p73m\nWEm6GJt37s17/txJzt6BXiQkgaxmMatpFi6fczmeOviUMSd+qC8NxRKcBDBpxs7R+wweWWFMrcxs\nmmnMx7P+vPwY2GnCchuBf23rSwAKg3uQihjTh8Xs3XcBksaWM2fjZCL/S7IIIPUn8FczkjjjnckF\n5YqNJw5jx54erFrSivaz23HH1hTey3mGU28Vvi/TB1AxVSb2+x06vmS0aRgAtLrU4q9bui4vxw4A\nyCZx6bTrQ92XSsdUDE04dvrFFKQBYFjVuJo0SG266cNiFt4B4NixyMHeqs7pkJ6Rl+oJsnI1ikVQ\nQPha/Paz23HV7K9CMy1QBbKDLUj3rsaDT87wTHFR9Dhip4oIUxES9X2d5Yomdu4895mC1KabAu8h\nnY458g5mDg2jN+n+n5wz8J/QBtw1tCbvmm7fGmylLoLKVUwt/q7nW/Fe3/q8Y0Mo/LOk8cXATmUX\nNrcdJdOI2k1PX3o0DWJbtaTV8xlNgfeuoTW4s+E+rDvWN5JjN39RbhiajIPZ6Zgt/4VDegbuGlqD\nR7OXYerkJC7d9AQO9aXRnEoiWSfIDBeuP6kTwe9e6P18Yfm9X6fx7JFDwTGwU9l5VZeMd2APG2By\nP3BMi2+A/JWXyy+6Hg8+OSPv/f207uO4ful8tL/+V8A7R7Fx+jT0i+RtAa/ZJPrevga/2Zg7Uflz\nTPqvJrx3fCmOnbBKLPvSGSQTgqmTkzh2IpPXoHdYFT/c3YO2s6YF/nOM+ptT0JW5NL648pTKbv56\n8244AuCNTd4TfKWubHTLrZu7m1taW1K4ZU26YGKwXuohIshkx+raG+sacdXsr2LX862ewfKcO76B\nhpxKk1NHVlivn7Udkhi7HrJJpHtX501e2s8EwPheWltSeHr9J7z+GACMzyraSq3MnSjKuvJURD4J\nYAuAOgD3quqmKK5L8bNjTw8EQN2UPXkldKeOrMCZiUs8XxvFysYbV5xrDDy/e2ErfvDsm8bXHOpL\nGxffDOlQwafByeGT2HbgLsz6wCx8Ydr12PV8K2546EVs3rl3NMDv2NOD7LtLMeAI1k3v35Qf1AEg\nkcGkGTsLArvXNw/T70wj845HX4n8m1PUPXKoOCUHdhGpA/BdAP8dwEEA/y4ij6rqL0u9NsXP5p17\nUTdlT97IVEY2fVh+1lmer3Vb2XjLz28BECy4J5tfxBm/8S30D76NbKYFkweuxq0f/6yVaunuHU15\n5JrcUBd6kU3vQC8efvdunMyuhmLJ6DxC1/6j+OHuHgwbvilL0tBy1+X4bI8R++yWVF4gb04lMTA4\nNJqX7+lL48ZtLxnz9EDp+fCweXmKXhQj9osA/FpV9wGAiDwIYCUABnYqcKgvjcnv31kwMpVEBk8f\n/T4A95pnt+Ca1WygkXveiH9kFyGkHkay+XwArQUrOW0Dg8OYk5xhdVQMQRyj7XRmGP/43AEMq6Le\n8I1FMy2QBkNwH8ovg8ytfDF9+7jivBl5x3NbINjcgjpQej68UhVPNCaKwN4K4EDOzwcBfDSC61IM\nzW5Jod9lZOo3KnZb2QgE2xPTb8Tfbxio2gG4f7Avb7ITADSbgIgA4l5l4xxt20Hd9I0l03chki27\nC3LsF075DH7dknINlM4gGqbyx8RZLhkmUFey4onGlG2BkoisFZEuEek6cqT4XVuoOrhtvuB3/vFZ\nf4KCCDnCb+m5X8+TwwOHPZ/Lb8Q/feYrecftAJxoyA/q9uKbk72fRvrQp6zNLlw4t7CrE7FG6oZv\nLPWnvYaTvauRHcxf4PP8ywtw44pz8camdjy9/hOjAdIZcK84bwY279zrWecexOade0cXFAXdCjD3\nteO11ysFF8WIvQfA3Jyf54wcy6Oq9wC4B7CqYiK4L1VI2ElM5/lil6DkxMMgO+/Y1775326BorCr\nYvOkZs/n8hvxN79vJ+Tw+aPzoaYADFjBeuD1nEU4sx8yXlMVo9Uutt/76Fzs6HPPpTuX7wPmBT6m\nkbHb5K9Jss7aeSmTLfxPMXeUHbY0tdh+OkzXRCuKEfu/AzhHRBaISAOA6wA8GsF1qUqFbc9qOh8C\nJCQReuuzTP+HMXT42oKOg0mZBFX1fC6/Ef/xzBF89uJ5o583QScz3TaWBgSNsx9C0/s3oX7KHrSk\nkvjGqoVoaXif8Wz36xQGxmLSLVMnJyGwyiE3f2oxNn96sWvfFzt4h11w5JafV1hlrktu34U/3/Fy\nqG8BFF7JI3ZVHRKRPwawE1a54/2q+orPy6iGhW3P6nZcVdH9+e5Q9968cy9O9C1G/XA2b/KxbuBq\nHJ/6A8/72x8ct/z8FmS1cMQ/s2kmvrFqIdrOmobNO/eiz2Uy0xmATx1ZgdSs7UDO6F4VELFGw3YO\nfXHTGQCAmy/+U9z0s7/I+zag2WTB6D6XM2AWU7myZ8PygmOrlrRiwfpOYw2/XVFjmnx1C+CmctJc\nx05kjN8syrVAbaKIJMeuqv+kqh9U1fer6jejuCZVr7C700e5m72dP65L7Yck+wEAkuzHcf1VoPu0\nn92OOy67w3MT5lVLWvH0+k/grmXrC85LyiTo0StRP2UPmt6/Caedtx6N79uJi6avGN18OSEJK92U\nQxIZ7D7+D9ixpwftZ7cj1X9dXi79pGERks3U/yVs5YrbyNzrWnaZpFMyIa79aFYtacXG1Qs97+cm\nzm0HwvT9jwK7O1JoQXentycyTXntYnaztxc3TTpzB5JTn4WIWm1vRdEw9VmcdfpZgZ6r/ex2dFzS\nMRqI3VJBpvP+12Vfx7UfmYvUyKSqiJWaef7oY+gd6MXMppnGbwMAgPq+0UnEWz/+WZx4fT3ee20T\nBl5fnxfUU8kEWltSo2kT06pNU+dFN4LCSpdcV5xXuOsSAGSGs8ayyNMa6/MmcJ0By/5gdJ9SNotr\n24GwE9BRYEsBKorf0n7nhGmuWU2zArcCyL0Phlpw4q3laJy9dTTNkSshCdxx2R3jvply0I2lnbKD\nLTjx+vrRtgkfuOMbBbXsQ8eXIDllDxZ88Cn0DhyGDLUg/dZyNGUuggjQdyIzOtkI5Jc6zj8jhWde\nPwoF8urkJyem47bL/qfrn4NXC2MTu/WDX/uAMNe1V/8++dqR2E2ouv05BG39kCtoSwEGdhpVah+W\n3NeLiHHk6ra7vdv1nB8Omk0CkilIddhe/vzLgZ+3WIu+twjq2lnGzP7PLDE8Fb+74Et4rLsX6eYH\nC/Lsmb4L0dCyOz9fn00WpGrcAiEAfPNnDxRcu7GuseBbiV2ZErY80g5IfgFrx54e3PjwS8bKm1x1\nIvi9j87FD3f3hOoxUyuVNW5zGEF6IxW8ppy9Yqj2ldqHxfl6twFDmKX5pmoaSWRcV4gmxJxZdAaA\n5Rf14Omj3y/6A8yrbNKN/UGk9cfw8P67oacnkTDUsienPgc4vo04V7AC1mTjA8++ORow7K/3G1cv\nxNQ5/4KTA/nXdi7gCtKXviWVxKmhbEGwtT9A/Cpm7CDb8egrxglY+3obVy8MXFaZ+2GU27itmhdC\nVaLjJXPsBKD0HeaNJY0GYSZMvT4ETMH90x/8dMExZ37zrewzeHj/3egd6IVCRz/A/BZY5QqzObSJ\nJDKQuhPm3xlSTIC59NJ5ph0Ig1QnBSmXFAGWzmtG3cinkgBICHDDQy/i0k1PoDmVNL4uN2CtWtKK\nF29bjv/c1I7/3NSOb1/7YeP8QZCyytz/L73ef7UJuxNVFDhiJwCl7zAf5LywE6ZuI2MZ/Z8xF8+8\nGH9+8Z8XnOsMYKZFR0HaEeTyK5scD1417rkO9aVxzgf8N5UOUoFy7EQGT79+dOwZYPXNAczNxwD/\ngOXWICzIqDbIh5HX+6pU6qYSHS85YicA3iWJQdoHuL1eciLwpLpJoZ7JdWRsyK/vf3c/gMJWB29n\nn8l/aYg+NZ37OnHR9y/Bwr9fiAv+fiEW3n8xvv7E9wFYwT1bwvyUDqdcU0oF5xpq3N0qTma3pAJV\nLY1HGmDq5GTRfdeDjGpL2e+1EpUpuexKIWdbiPHCwE4A3EsYL59zOTqe6fBNXZheXy/1qE+MfSns\nH+wPlfY5SzliAAAOJUlEQVRwlhu2THIftdp9YpzP2jhrO+qn7Bk9z23k6/xg6tzXiVv+7Vaks+8C\nMpIjrxvAw/v/N77+xPexY08PNNMc6H04aTaJU29dAx32Dq6mGnc7hfHZi+cZA+EV583AHVtT6Htz\nJWRoKuBSzhmmXDKoyQ31RQes3Pp3tzJPvw8jr28LE62HDatiaJSpKmbLC1uMX+tN1S3O16eH0ug7\nVThCzn1t0Eocr/JJ+5oAjM+qCs9dihrrGrHyAyvx1MGnfJ8dAGRoKk4/0oG3ss8U7njkIfc5ho4v\nKejyaDr/vdfG9qxxlseZmoA5K0uSdYKmhnr0pzMFKQDn6wdODblOcgbV6tGFMii3lIlpwteeQG31\nuV+UlSmVxHJHioRbaZ9AfNsBeJUFCgRTGqbgxNCJgq3l7NFlkPJJAEA2iU+ddQN+eGCzZxmiZpNI\n9V+HqxbNyquKuXzO5dj+q+15z+FFFRh4bVNBvbgOTwagSNQXpgySMgkne1fjxLHFecfrp+xxrcvP\nDo41GwuyvVyQunGv6wSplPHi3F6wmC3x/Grji82TR1lLXkksd6RIuE1gBqlu8SoLVCj6B/sLjudW\n4gQpn1QFTvauxoP7ZuCM3/DeDEMSGZyc+n08fTR/gdTHHvxY4KAOWLXo9mSfsxujHShM30Qy/R8u\nKP0bOr4E2bqE1Y9dT40eT8ok1A1cjRMAWiYnoYqCLfacguSgvXqymCb5rjhvBp587UigWne3KpUw\ngd2v7LHY3ZnctkQcz8qUSmJgJ0/rlq4rSIEErW4xvTaIwwOHA5dPaqbFCq4Yxqm3V6Bx+sO+r3PW\n6LulXKxGXo5j2Tp8asGXsHipd6BoP7vdmFJyjjqbU0nIqTa8ezCL1Jm7oPV9mNU0E5dOux67DrVC\nkUbfiUxevfYND72Irv1H8Y1VC/Ou7VZZ4uT1AeAVOD/89V2hUzVh+7+4nd/Tlx5tV1CMibYXKwM7\nebKDUzErUp2vDbpac2bTzEDlk9bk4hmjP79z+Hz89TXnu84L5Apao58dSkHqRuqmhyfj1FtX47bf\nH9u+r5hAYQfP/LTDEmSOL0EqWYeLLmzFg0/2IJ0x12srgAeefRNtZ03Lu59fZ0VbsRUxHdecHzpV\nE/ZeXh9OpS5Amkh7sTLHTmVz2T9eZky/5LJz7EGCMwCoCt57bSOA/Hyp32QrYOX5vdJF2aEUBn51\nW96xKHOybnnfOhHjZtdOpmfJ/TbQMjmJ904O5S3pLybvbbp+kG8GYe+1Y08Pvv7jV4wbittqLSce\nNebYqap07uvEiSHzastcKz+wcnSkv/7f1vucDQAje4i+byeOJ/uxfFv+Nwqva3h9M1AFTr11Td6x\nZJ17u9piuKUdggR1t9c7R6VRL8qxr+/1oZRVDX2voBO3cW7tGyUGdiqLLS9sCTRB+dTBpwBYaZxg\ngR15m1w48+deI3+vcs6m+imoz1yEPljXnTo5iduuPj/Sr/JBc+Jer/czXukHt8nI3BG63dI3yIdK\n0B2h4traN2oM7FQWxbQmmNU0yzcdIwJAClsEbHxuI9rPbnedwL323GtHR/WmyeENl9yC9s8V7jhU\nDLdRsyk4OksG3QSt6BivZfR+k5GmPVm9cuRBRuJxrmKJGgM7lUXQjoi5ZZTFVtUA1irXzn2dvpO/\nvpPD3VuBx28H+g8CzXOAZRuARWvy7uW1yCpIgMsNjl4j+LCLf8IG1zD8PjDcyhY7Hn3F+Dq3915s\namei4+QpjYqyH7vz9UEmM509wzv3dWLjcxt9J1zdhOn9btS9FfjxV4FMTsBJpoCrvzMa3E3vK/d9\nhF0YE+VCmqgX5bi1zAUK0zBuKz2d7NcB8E3tUPDJU/aKIQAw9lkJ09fF7/WmbeauPfda1+3p7Ovl\nBvV6CfcFM0zvd6PHb88P6oD18+O3j/7o1+7Yqy7btAdmlC1eg7TCDSpsy9ygufDcxUd+vWIoOKZi\nIlbqqLdS3ALUpuc3BXp+rwCXm/YI+mdhut6QDiEhicCtcovZLDtP/0Hf437tjt1SDIKx1rd2iqRr\n/1E8+doRpDPDoyWPfj1QvES5wUPYlrlBa+pzXzeR6szHG0fsESp11Ot1Xb+2uaVyC1B9p/oC3S9o\nP/eg78XtelnNGrtIJhP5mz4Us1l2geY5vse92h0D5hG4aYLU3hHJDsTDqqMj9WKDXTlG/7mcG2w4\nR+BTJ/tvzEHR4Ig9QkFGrWEF3bKu1G8KXpObzue379U70Ds6gk5IwtjPZUrDFCzfthyHBw4XNP3y\n2n7P7XnsjbCd79V+zki/KS3bYM6xL9sw+qNfy4UwE6RR9FrJFeUyer+JXdMHhqmmfiL1a6mkkiZP\nRWQzgKsBDAJ4HcDvq6q58UaOuE6eltIJ0c3ybct92+b6TeAF0bmv07VuPPf5g0yC2uqlHiLiW7/u\n1gK41PcUiRKrYkyCdGG0VUtb2VJa5jqvM1H6tYyHsrTtFZHlAJ5Q1SERuRMAVPUmv9fFNbAHCcJh\nBfmwiOq+bkv+E5LAHZfdgfaz213vlXuuqvr2NHd7L7lqdb7Cj1eQdKqmJfQMypVXlpYCqpobNZ4F\n8KlSrlfrSumE6CZI21y3fHTvQC8WfW9R4KB480dvNo7Gs5odTZn4VZqo6miQXvS9RZ7n2ky7F+UG\n9I0f2xiLgG5za4/r3CSj2tIUnNysHVHm2L8I4CG3X4rIWgBrAWDevHkR3rZ6lNIJ0U2QDwu/vude\nuWzT85s2abbnCvwWGuUG6SCLkpzvJeicQq0zBcm2s6ZxREyR8E3FiMi/ADBN/d+qqo+MnHMrgDYA\nqzVAbieuqZjx4peSCJr3Dpqa8Ur/bPzYRtd7mRYY+T3Xpo/ll1OORzqLKC4iS8Wo6m/53OgLAK4C\nsCxIUKfw/Oq/g/Y9D7pgxyv9k3uv3KoYu1ol9zm9vgEAVrB2vq+gZZNE5K6kVIyIfBLAnwH4uKr6\n92SlsnBbxBN0wY5f+ifMQiOvRlumuYdStuIjIkupC5T+GsDpAH4qIi+KyN9G8EwUknNhlCmoh5nE\nNS3/L6XMMMz11i1dV7AAKZLFRkQTCJuAxYBbXjq39LCWSgXjWuZIVCruoFRjSglm7rsAadELo4pV\nzPswvYYTpUTFY2CvAqWW+FVLXrqY9zFRyhuJyolNwKqAX+tXP9WSly7mfZT63omoEEfsVaDUEr/x\nWBhVjGLeB8sbiaLHwF4FokilhClBHC/FvI9qSSMRxQlTMVWgWlIppSrmfcTlvRNVE47YixB1Od54\np1LKVT5YzPuoljQSUZywjj2kqukTHlCtPS8RuYvtZtbl2CbOS61VcdTa8xJR6WoqFVMNNc+1VsVR\na89LRKWrqRF7NYw+/TYvrja19rxEVLqaCuzVMPqstSqOWnteIipdTaViqqHmudaqOGrteYmodDVV\nFcMKDyKayGLZ3ZGjTyIifzUV2IHqWDpPRFTNamrylIiI/DGwExHFDAM7EVHMMLATEcUMAzsRUcww\nsBMRxQwDOxFRzDCwExHFTEVaCojIEQD7y3zb6QDeKfM9o8Tnr6xafv5afnaAz5/rLFWd4XdSRQJ7\nJYhIV5AeC9WKz19Ztfz8tfzsAJ+/GEzFEBHFDAM7EVHMTKTAfk+lH6BEfP7KquXnr+VnB/j8oU2Y\nHDsR0UQxkUbsREQTwoQL7CLyFRF5TUReEZG7Kv08xRCRr4mIisj0Sj9LGCKyeeTPvltEfiQiLZV+\nJj8i8kkR2SsivxaR9ZV+njBEZK6IPCkivxz5+15zG92KSJ2I7BGRxyr9LGGJSIuIbBv5O/+qiPy3\nct17QgV2EbkCwEoAi1X1fAB/WeFHCk1E5gJYDuDNSj9LEX4K4AJVXQTgPwDcXOHn8SQidQC+C+BK\nAB8C8Hsi8qHKPlUoQwC+pqofAnAxgC/X2PMDwDoAr1b6IYq0BcA/q+p5ABajjO9jQgV2AH8IYJOq\nngIAVX27ws9TjLsB/BmAmpscUdVdqjo08uOzAOZU8nkCuAjAr1V1n6oOAngQ1sCgJqhqr6q+MPLv\n78IKLK2VfargRGQOgHYA91b6WcISkWYAlwO4DwBUdVBV+8p1/4kW2D8I4GMi8pyI/ExEPlLpBwpD\nRFYC6FHVlyr9LBH4IoCfVPohfLQCOJDz80HUUGDMJSLzASwB8FxlnySUb8MaxGQr/SBFWADgCIC/\nG0kl3SsiTeW6ec3teepHRP4FwEzDr26F9X6nwfpa+hEAW0XkbK2i0iCf578FVhqmank9v6o+MnLO\nrbDSBA+U89kmKhE5DcAPAfyJqh6v9PMEISJXAXhbVXeLyG9W+nmKUA9gKYCvqOpzIrIFwHoAf1Gu\nm8eKqv6W2+9E5A8BbB8J5M+LSBZWH4cj5Xo+P27PLyILYY0CXhIRwEpjvCAiF6nq4TI+oievP38A\nEJEvALgKwLJq+kB10QNgbs7Pc0aO1QwRScIK6g+o6vZKP08IlwK4RkR+G0AjgCki8gNV/VyFnyuo\ngwAOqqr9DWkbrMBeFhMtFbMDwBUAICIfBNCAGmkupKovq+r7VHW+qs6H9RdnaTUFdT8i8klYX62v\nUdUTlX6eAP4dwDkiskBEGgBcB+DRCj9TYGKNAO4D8KqqfqvSzxOGqt6sqnNG/q5fB+CJGgrqGPnv\n8oCInDtyaBmAX5br/rEbsfu4H8D9IvILAIMAPl8Do8Y4+WsAkwD8dORbx7Oq+geVfSR3qjokIn8M\nYCeAOgD3q+orFX6sMC4FcD2Al0XkxZFjt6jqP1XwmSaSrwB4YGRQsA/A75frxlx5SkQUMxMtFUNE\nFHsM7EREMcPATkQUMwzsREQxw8BORBQzDOxERDHDwE5EFDMM7EREMfP/AZ05j0R4u4OkAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3b0d97ba58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.utils import shuffle\n",
    "\n",
    "n_samples_per_component = 100\n",
    "\n",
    "samples = []\n",
    "component_ids = []\n",
    "\n",
    "for i, mean, Cov in zip(range(n_components), true_means, true_covariances):\n",
    "    data = rng.multivariate_normal(mean, Cov, size=n_samples_per_component)\n",
    "    samples.append(data)\n",
    "    component_ids.append(i * np.ones(n_samples_per_component, dtype=np.int32))\n",
    "    plt.scatter(data[:, 0], data[:, 1])\n",
    "    \n",
    "data = np.vstack(samples)\n",
    "component_ids = np.concatenate(component_ids)\n",
    "\n",
    "data, component_ids = shuffle(data, component_ids, random_state=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no closed form formula for the MLE. Let's use the scikit-learn implementation of the EM algorithm instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "\n",
    "gmm_em = GaussianMixture(n_components=3, random_state=0).fit(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Average loglikelihood of the data under the EM-MLE model (higher is better):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3.6213242571063557"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm_em.score(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.60843729, -0.7027044 ],\n",
       "       [ 1.94313479,  4.48270848],\n",
       "       [ 1.51596953, -0.53611044]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm_em.means_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([ 1.49014246, -0.4147929 ]),\n",
       " array([ 1.94306561,  4.56908957]),\n",
       " array([-0.70246012, -0.70241087])]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "true_means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.32384082,  0.35908973,  0.31706944])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm_em.weights_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 3.7638645 ,  1.71761884],\n",
       "        [ 1.71761884,  1.15973123]],\n",
       "\n",
       "       [[ 0.30253199,  0.76447466],\n",
       "        [ 0.76447466,  4.88339622]],\n",
       "\n",
       "       [[ 2.72596464, -0.17407796],\n",
       "        [-0.17407796,  0.37998651]]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm_em.covariances_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 3.08286918, -0.32502055],\n",
       "        [-0.32502055,  0.5147776 ]],\n",
       "\n",
       "       [[ 0.43166016,  0.77894194],\n",
       "        [ 0.77894194,  3.71918704]],\n",
       "\n",
       "       [[ 3.29150879,  1.5703531 ],\n",
       "        [ 1.5703531 ,  1.12457826]]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.asarray(true_covariances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's find the MLE by gradient descent. First we need a helper function to compute the log of the sum of likelihoods of the components in a numerically stable way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " 0.6123\n",
       " 1.8312\n",
       " 1.1842\n",
       " 0.6889\n",
       "[torch.FloatTensor of size 4]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def logsumexp(data, dim=0):\n",
    "    \"\"\"Numerically stable log sum exp\"\"\"\n",
    "    max_score, _ = torch.max(data, dim)\n",
    "    if dim == 0:\n",
    "        max_score_bcast = max_score\n",
    "    elif dim == 1:\n",
    "        max_score_bcast = max_score.view(-1, 1)\n",
    "    else:\n",
    "        raise NotImplemented(\"logsumexp with dim=%d is not supported\" % dim)\n",
    "    return max_score + torch.log(torch.sum(torch.exp(data - max_score_bcast), dim))\n",
    "\n",
    "\n",
    "test_data = torch.randn(3, 4)\n",
    "logsumexp(test_data, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " 0.6123\n",
       " 1.8312\n",
       " 1.1842\n",
       " 0.6889\n",
       "[torch.FloatTensor of size 4]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "torch.log(torch.sum(torch.exp(test_data), 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = Variable(torch.FloatTensor(data))\n",
    "\n",
    "\n",
    "logsotfmax = torch.nn.LogSoftmax()\n",
    "weights = Variable(torch.ones(1, n_components), requires_grad=True)\n",
    "\n",
    "means = []\n",
    "prec_diags = []\n",
    "prec_off_diags = []\n",
    "for i in range(n_components):\n",
    "    mu = Variable(torch.randn(n_features), requires_grad=True)\n",
    "    means.append(mu)\n",
    "    d = Variable(torch.ones(n_features), requires_grad=True)\n",
    "    prec_diags.append(d)\n",
    "    W = Variable(torch.randn(n_features, n_features), requires_grad=True)\n",
    "    prec_off_diags.append(W)\n",
    "\n",
    "\n",
    "def mixture_nll(X, weights, means, prec_diags, prec_off_diags):\n",
    "    log_normalized_weights = logsotfmax(weights).transpose(1, 0)\n",
    "    logliks = []\n",
    "    for mu, d, W in zip(means, prec_diags, prec_off_diags):\n",
    "        logliks.append(loglik(X, mu, d, W))\n",
    "    \n",
    "    logliks = torch.cat(logliks).view(n_components, -1)\n",
    "    return torch.mean(-logsumexp(logliks + log_normalized_weights, dim=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 42.4368\n",
       "[torch.FloatTensor of size 1]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mixture_nll(X, weights, means, prec_diags, prec_off_diags)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params = [weights]\n",
    "params.extend(means)\n",
    "params.extend(prec_diags)\n",
    "params.extend(prec_off_diags)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       " 1  1  1\n",
       "[torch.FloatTensor of size 1x3]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 42.43677520751953\n",
      "100 3.063232660293579\n",
      "200 2.824497938156128\n",
      "300 2.7710728645324707\n",
      "400 2.750455141067505\n",
      "500 2.7379977703094482\n",
      "600 2.7285289764404297\n",
      "700 2.720960855484009\n",
      "800 2.7149853706359863\n",
      "900 2.7103676795959473\n",
      "1000 2.706846237182617\n",
      "1100 2.7041540145874023\n",
      "1200 2.7020885944366455\n",
      "1300 2.7006301879882812\n",
      "1400 2.6996991634368896\n",
      "1500 2.699103593826294\n",
      "1600 2.6987264156341553\n",
      "1700 2.698498487472534\n",
      "1800 2.6983683109283447\n",
      "1900 2.698298931121826\n",
      "1982 2.6982688903808594 converged!\n"
     ]
    }
   ],
   "source": [
    "learning_rate = 0.05\n",
    "optimizer = torch.optim.Adam(params, lr=learning_rate)\n",
    "best_loss = np.inf\n",
    "for t in range(5000):\n",
    "    loss = mixture_nll(X, weights, means, prec_diags, prec_off_diags)\n",
    "    if t % 100 == 0:\n",
    "        print(t, loss.data[0])\n",
    "    if loss.data.numpy() < best_loss:\n",
    "        best_loss = loss.data.numpy()\n",
    "    else:\n",
    "        print(t, loss.data[0], 'converged!')\n",
    "        break\n",
    "    optimizer.zero_grad()\n",
    "    loss.backward()\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([ 1.94653308,  4.54145813], dtype=float32),\n",
       " array([-0.34557292, -0.51862174], dtype=float32),\n",
       " array([ 2.11376023, -0.76199687], dtype=float32)]"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[m.data.numpy() for m in means]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.60843729, -0.7027044 ],\n",
       "       [ 1.94313479,  4.48270848],\n",
       "       [ 1.51596953, -0.53611044]])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm_em.means_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([ 1.49014246, -0.4147929 ]),\n",
       " array([ 1.94306561,  4.56908957]),\n",
       " array([-0.70246012, -0.70241087])]"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "true_means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.3512145 ,  0.43690208,  0.21188341], dtype=float32)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "softmax = torch.nn.Softmax()\n",
    "softmax(weights).view(-1).data.numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.32384082,  0.35908973,  0.31706944])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gmm_em.weights_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
